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Abstract 

We perform some experimental simulations in spherical symmetry and axisymmetry to understand 
the post-shock-revival evolution of core-collapse supernovae. Assuming that the stalled shock wave is 
relaunched by neutrino heating and employing the so-called light bulb approximation, we induce shock 
revival by raising the neutrino luminosity by hand up to the critical value, which is also determined 
by dynamical simulations. We incorporate nuclear network calculations with a consistent equation of 
state in the simulations to account for the energy release by nuclear reactions and their feedback to 
hydrodynamics. Varying the shock-relaunch time rather arbitrarily, we investigate the ensuing long- 
term evolutions systematically, paying particular attention to the explosion energy and nucleosynthctic 
yields as a function of this relaunch time, or equivalently the accretion rate at shock revival. We study 
in detail how the diagnostic explosion energy approaches the asymptotic value and which physical 
processes contribute to the explosion energy in what proportions as well as their dependence on the 
relaunch time and the dimension of dynamics. We find that the contribution of nuclear reactions 
to the explosion energy is comparable to or greater than that of neutrino heating. In particular, 
recombinations are dominant over burnings in the contributions of nuclear reactions. Interestingly 
ID models studied in this paper cannot produce the appropriate explosion energy and nickel mass 
simultaneously, overproducing nickels, whereas this problem is resolved in 2D models. 
Subject headings: Nuclear reactions, nucleosynthesis, abundances — stars: supernovae: general — 
MHD — methods: numerical 



1. INTRODUCTION 

The mechanism of core-collapse supernovae (CCSNc) 
has refused our understanding for the last half century. 
CCSNe are initiated as the implosion of central cores of 
massive stars at the end of their lives. The inward motion 
is halted when the density exceeds the nuclear saturation 
density, ~ 3 x 10 14 g/cm 3 , above which matter becomes 
very stiff with the adiabatic index of Y > 2. At this 
point, the core is still lepton-rich with the lepton frac- 
tion of Yi ~ 0.3 and very hot with the temperature of a 
few tens MeV compared with the ordinary neutron star, 
which is the end product of CCSNe. The gravitational 
energy liberated so far is ~ 10 53 erg, far greater than the 
typical kinetic energy of ejecta in supernova explosions, 
~ 10 51 erg. This large energy is mostly stored in the core 
as internal energy and unavailable for the prompt explo- 
sion initiated by the shock wave produced at the core 
bounce. As a result, the shock wave is soon stagnated 
in the core owing to dissociations of nuclei and neutrino 
cooling. The current study on the mechanism of CCSNe 
is hence focused on the revival of the stalled shock wave. 

There ar e a couple of viab l e mechanism s proposed 
at present (jKotake et alj|2012t Uankall2012[ ). The most 
promising among them is supposed to be the neutrino 
heating mechanism. The large energy reservoir in the 
core is tapped by neutrinos over the next ten seconds. 

1 Science & Engineering, Waseda University, 3-4-1 Okubo, Shin- 
juku, Tokyo 169-8555, Japan 

2 Kumamoto National College of Technology, 2659-2 Suya, 
Goshi, Kumamoto 861-1102, Japan 

3 Yukawa Institute for Theoretical Physics, Kyoto University, 
Oiwakc-cho, Kitashirakawa, Sakyo-ku, Kyoto, 606-8502, Japan 

4 Advanced Research Institute for Science & Engineering, 
Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan 



In this scenario, a part of the energy carried away by 
neutrinos are deposited in the so-called heating region, 
where the heating of matter via absorptions of neutri- 
nos emitted by the proto neutron star overwhelms the 
cooling through local emissions of neutrinos. If the en- 
ergy deposition is large enough, the accretion shock wave 
will be reinvigorated and restart to propagate outward, 
expelling the stellar envelope and culminating in well- 
known optical displays of supernovae when the shock 
wave reaches the stellar surface. In fact, it is now al- 
most established that there is a critical luminosity for 
the revival of the stalled shock wave, which is a function 
of the rate of matter accretion: the more rapidly mat- 
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existence of the critical luminosity, calculating a sequence 
of spherically symmetric, steady, shocked accretion flows 
and showing the non-existence of such flows for luminosi- 
ties beyond a certain threshold giv e n for each mass ac- 
cretion rate. lYamasaki fc Yamadal (|2005[ ) analyzed the 
sequence more in detail and showed that there are in 
general two accretion flows for a given accretion rate, 
with one of them with a larger shock radius being al- 
ways unstable to radial perturbations and the other be- 
ing stable and that the br anches join with each othe r 
at the c ritical luminosity (iPejcha fc Thompson! |2012| ) . 
Recently iPeicha fc Thompson! (j2012t ) found that the so- 
called antcsonic condition predicts accurately the point, 
at which the spherically sy mmetric, steady, sho cked ac- 
cretion flow ceases to exist. lOhnishi et al.1 (|2006l ) demon- 
strated for the spherical case that the accretion shock 
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may be revived by overstabilizing radial oscillations be- 
fore the point of the non-existence of steady sh ocked flows 
is reac hed, the fact also confirmed recently by IFernandezI 
(|20ll . 

For more realistic, non- spherical flows, the criti- 
cal l u minosity is lowered (1 Jankal 119961: lOhnishi et alJ 
20061: IMurphv fe Burrowi 120081: INordhaus et alJ |2010T 



Hanke et al.l I2012TT Non-radial instabilities referred to 



as the standing accretion shock instabilities or SASI are 
known to set in for lower neutrino luminosities than the 
radia l instabilities (Foglizzo 2001J lYamasaki fe Yamadal 
|2005[ ). They effectively enhance the neutrino heating 
by dredging up heated material and channeling cold 
material. It has been demonstrated that the ratio of 
the residence time to the heating time is a good mea- 
sure to gauge ho w close to sh o ck revival a par t icular 
accretion flow is dJanka l2001t iThompson et al.l 120031 : 
IMurphv fe Burrowsl 12008 ) . Here the residence time is 
a duration, in which each mass element lingers in the 
heating region and the heating time is a time needed 
to deposit enough energy to unbound the mass element. 
This criterion is rather easy to evaluate and has been ap- 
plied particularly for numerical r esults dMarek fe Jankal 
20091: IMurphv fe Burrowsl l2008t INordhaus et alJ 120101: 
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It should be emphasized that the revival of stalled 
shock waves is not sufficient to ascertain the supernova 
mechanism. In fact, almost all numerical simulations 
that have reported successful shock revival have yielded 
energies of ejecta that are substantially smaller than the 
canonical value of 10 51 crg. Although the authors claim 
that the energies are still increasing at the end of their 
simulations and might reach appropriate values should 
the computations be continued long enough, a convinc- 
ing demonstration remains to be done. In addition to the 
explosion energy, the mass of neutron star as well as nu- 
clcosynthetic yields should be also reproduced properly 
in the successful supernova simulations. 

In this paper, we are interested in what will happen 
after the stagnated shock wave is successfully relaunched 
by the neutrino heating. In particular we discuss (1) 
when the explosion energy is determined, (2) which pro- 
cesses contribute to the explosion in what proportions, 
(3) how the explosion energy and neutron star mass are 
dependent on the timing of shock revival, and (4) how 
multi-dimensionality affects all of these issues. For these 
purposes, we have done a couple of numerical experi- 
ments in ID (spherical symmetry) and 2D (axisymme- 
try). Controlling neutrino luminosities under the light 
bulb approximation, we have induced shock revival from 
different points on the critical curve (the critical lumi- 
nosity as a function of the mass accretion rate) and com- 
puted the following evolutions of matter flows outside the 
proto neutron star long enough for the energy of ejecta 
to become constant. In so doing we have taken into 
account nuclear reactions in a manner consistent with 
the EOS employed. The feedback from the reactions 
to hydrodynamics is fully incorporated . These consis- 
tencies were lacked in previous works dMarek 
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The paper is organized as follows. In the next section, 
we describe the models and numerical details together 
with the assumptions and approximations employed in 



this paper. The main results are presented in Section [3] 
first for the spherically symmetric ID case and then for 
the axisymmetric 2D case. The summary and some dis- 
cussions are given in the last section. 

2. SETUP 

2.1. Outline 

Before going into the details of our model building, 
we give a brief description of what we are going to do, 
emphasizing the underlying ideas. 

We are interested in what happens after the relaunch 
of the stalled shock wave. The investigations in this pa- 
per are of experimental nature. We assume that the 
neutrino heating mechanism works successfully, which 
implies that the neutrino luminosity and accretion rate 
should be located on the critical curve just at the shock 
revival. Exactly at which point on the curve the shock is 
relaunched is still uncertain as mentioned earlier, how- 
ever. We hence take the neutrino luminosity (or equiv- 
alently the accretion rate) at the shock revival as a free 
parameter and vary it arbitrarily to see how the ensu- 
ing physical processes are affected. We prepare a cou- 
ple of initial conditions, which correspond to the points 
for different neutrino luminosities on the critical curve. 
We then solve the hydrodynamical equations together 
with nuclear reactions in ID (spherical symmetry) and 
2D (axial symmetry) to obtain the ensuing evolutions. 
See Section 12.41 for more details on how to trigger the 
shock revival. 

We do not solve the evolution of the central high- 
density region, in which a proto-neutron star sits, but 
replace it by appropriate inner boundary conditions. The 
temporal variation of mass accretion rate is obtained by 
the computation of the infall of a realistic stellar envelope 
after the loss of pressure support at the inner boundary 
and is employed for the preparation of the initial states 
and subsequent hydrodynamical simulations. This also 
enables us to use the mass accretion rate as a clock. We 
follow the post-revival evolutions long enough so that 
the so-called diagnostic explosion energy is settled to 
the terminal value. Integrating the heating rates both 
by neutrino absorptions and nuclear reactions, we also 
obtain each contribution to the explosion energy. The 
nuclear reactions are divided into the recombinations of 
free nucleons to heavy nuclei and the explosive nuclear 
burnings and we estimate their contributions separately. 
We further distinguish the recombinations that occur in 
the nuclear statistical equilibrium (NSE) and those out 
of equilibrium. By so doing, we can pin down what con- 
tributes to the explosion energy in what proportions. We 
can also find the dependence of the results on the tim- 
ing of shock relaunch as well as the dimensionality of 
dynamics. 

In the following we give the details of our modeling. 
The basic equations, input physics and numerical meth- 
ods are described first. Then the preparation of the ini- 
tial conditions, which requires multiple steps to avoid full 
computations of the collapse to bounce to shock stagna- 
tions, will then be presented in detail. 

2.2. Basic equations, input physics & numerical 
methods 

Throughout this paper we neglect relativity and em- 
ploy Newtonian equations of motion. Not only for the 
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post-relaunch evolutions but also for the preparations of 
the initial conditions we solve the following equations: 
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where p, P, v, e, Y e , Y{ and Na are mass density, pres- 
sure fluid velocity, energy density, electron fraction, num- 
ber fraction of nucleus i and Avogadro's number, respec- 
tively. We denote the Lagrange derivative as D/Dt. Note 
that the energy density in Eq. ((3]) includes the rest mass 
energy and the energy production by nuclear reactions 
are thus taken account. 

In Eq. ((2]), it is expressed explicitly that the gravi- 
tational potential has two contributions, $ from the ac- 
creting matter and $ c from a central object, whose mass, 
Mj n , is a function of time and calculated by the integra- 
tion of mass accretion rates at the inner boundary of 
computational domain. They satisfy the following equa- 
tions: 

A $ = AirGp, (6) 

and 



GMh 



(7) 



where G is the gravitational constant. 

Q v and N v in Eqs. ((3|) and (HJ) are the source terms 
that give the rates of the changes in specific energy den- 
sity and electron fraction, respectively, owing to weak 
interactions. In the present study, we take into account 
only absorptions and emissions of electron- type neutrinos 
and anti-neutrinos on nucleons. The wea k interaction 
rates are adopted from IScheck et al.l (|2006|) . We employ 
in this paper the so-called light bulb approximation, in 
which the computation of neutrino transport is neglected 
and neutrinos are just assumed to be emitted isotropi- 
cally from the neutrino sphere with a given luminosity 
(Eq. (HU])) and energy spectrum, which we assume t o 
have the Fermi-Dirac distribution (|Ohnishi et al.ll2~006l ). 
The radius of neutrino sphere is a function of time and is 
assum ed in this paper to be given as follows (jScheck et al.l 
2006): 



1 + (1 - exg{-t/-to)){Ru,i/Ru,f - 1) 



(8) 



where and R v j are the initial and final values, re- 
spectively, and t c is the characteristic time scale. They 
are set to be R v ^ = 58km for v ei 52km for v ei R v j = 
15km and t c = 800ms. The temperatures, T„, in the 
Fermi-Dirac distributions for the electron-type neutrino 
and anti-neutrino are chosen so that their average ener- 
gies would be (e^) = 20 - 8.0 x (l/2) t / 200ms MeV and 
(e, 7 „) = 23 - 8.0" x (l/2)*/ 200ms MeV dSumivoshi et al.l 
l2005i ). The chemical potentials are assumed to be zero. 



In evaluating the heating and cooling of matter by neu- 
trino absorptions and emissions, we employ the local dis- 
tribution function of neutrino given by the following ex- 
pression: 

f(r,e)= C{ f) (9) 

1 + cxp(e/fc B r„) 

where ks is the Boltzmann's constant and the normaliza- 
tion factor, C(r), is determined so that the local number 
density of neutrino, n v (r), is given by the following rela- 
tion: 

L„=4Trr 2 n„{r)-(e v )-(lM(r)), (10) 

where the last factor, (p(r)), is the flux factor that ac- 
counts for how quickly the angular distribution becomes 
forward -peaked. We again employ the fitting formula 
given in lScheck et a l. (2006) for the radial dependence of 
the flux factor. 

Eq. ([5]) expresses the nuclear reactions. We deploy 
28 nuclei: n, p, D, T, 3 He. 4 Hc and 12 a-nuclci, i.e., 
12 C, 16 0, 20 Ne, 24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 
52 Fe, 56 Ni, and 9 their neutron-rich neighbors, or 27 A1, 
31 P, 35 C1, 39 K, 43 Sc, 47 V, 51 Mn, 53 Fc, and 55 Co. Wc 
take into account emissions of a nucleon and a particle 
as one-body interactions as well as three-body reactions 
such as 3a — > C in addition to the main reactions: (a, 7), 
(a,p), (p, 7) and their i nverses. The reaction rates ar e 
taken from REACLIB (jRauscher fc Thielemann] l2000f ). 
As demonstrated later, the employment of this rather 
small network is validated by the re-computations of nu- 
clear yields with a larger network including 463 nuclei , 
from n, p, and 4 He up to 94 Kr (|Fuiimoto et al.ll2004h 
for the densities and temperatures obtained by the simu- 
lations as a post-pro cess. We comb ine the semi-implicit 
scheme developed bv lTimmesI ()1999l ) and the full implicit 
method to numerically handle Eq. ([5]). The network cal- 
culations are needed only at low temperatures, since the 
NSE is achieved at high temperatures and the composi- 
tions are determined so that the free energy should be 
minimized. In fact, the network computations are ex- 
pensive when the NSE is established. We hence solve 
the network only for T < 7 x 10 9 K and otherwise calcu- 
late the NSE compositions for the same 28 nuclei as used 
for the network computations. This temperature is high 
enough to ensure the establishment of NSE. 

The equation of state we use in this paper is the sum 
of the contributions from nucleons, nuclei, photons, elec- 
trons and positrons. The first two are treated as ideal 
Boltzmann gases, the composition of which is obtained 
either by the network computations or by the NSE cal- 
culations as mentioned above. The photons are an ideal 
Bose gas and easy to handle whereas the electrons and 
positrons are treated as ideal Fermi gases, in which ar- 
bitrary degeneracy and relativistic kinem atics are fully 
taken into account (|Blinnikov et al.l lT996). It should be 
repeated that we include the rest mass contribution in 
the energy density so that the energy release by nuclear 
reactions is automatically taken into account properly. 

The hydrodyn amical evolutions are s o lved with the 
ZEU S-2D code (|Stone fc Normanl 119921 : lOhnishi et all 
|2006| ) on the spherical coordinates. For the ID compu- 
tations, we use the same code just suppressing motions 
in the 9 direc tion. We adopt a 15MQ progenitor model 
computed by IWooslev fc Hegerl (|2007l ) . The profile just 
prior to collapse is shown in Fig. [TJ In this paper we 
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Fig. 1. — The initial profile of the 15Mq progenitor. The density 
and temperature as well as mass fractions of some representative 
nuclei are shown. 

concentrate on this canonical model. The dependence 
of the results on the progenitor structures and the de- 
tails of our modeling will be reported els ewhere. Note 
that a ccording to the recent observations (jSmartt et al.l 
I2009al ISmartdl2009rj ISmith et al.ll2011[ ) of core-collapse 
supernovae as well as their progenitors, the 15M Q star 
may be a typical progenitor of Type-II supernovae. 

2.3. Step 1: ID simulation of the injall of envelope 

We now proceed to the description of the first step 
in the preparation of initial models. The aim of this 
step is to sample the mass accretion rates as a function 
of time together with the changes in the structure and 
composition of the envelope. For this purpose we perform 
a ID simulation of the spherically symmetric implosion 
of the stellar envelope. We excise the interior of r = 
60km and replace it with the inner boundary, at which 
we impose the free inflow condition. Since no bounce 
occurs in this simulation, no shock wave emerges. Note 
that what we need is the accretion rate and the structure 
outside the stalled shock wave, which would be produced 
and stalled somewhere inside the core in reality and that 
they arc unaffected by what happens inside the shock 
wave, since they are causally disconnected. The location 
of the inner boundary is chosen so that they would always 
reside inside the shock wave. 

We deploy 500 grid points to cover the region extend- 
ing up to r = 2 x 10 5 km. This is large enough to ensure 
that matter outside the outer boundary does not move 
essentially for ~ a second during this stage. The weak 
interactions are turned off for this computation, since 
the they are indeed negligible in the infalling envelope. 
Note again that the computational results for the region 
that would be engulfed by the shock wave in reality are 
irrelevant and do not have any consequence on the re- 
sults outside. Hence the neglect of neutrino heating and 
cooling is completely justified. The nuclear reactions for 
the 28 nuclei, on the other hand, are computed for the 
region with T < 7 x 10 9 K to follow the change in chemi- 
cal composition and take account of its influence on the 
hydrodynamics during the implosion. The NSE compo- 
sition is calculated for higher temperatures. 

We show the results in Fig. [2] In the upper panel we 
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Fig. 2. — The mass accretion rates. The upper panel shows them 
as a function of radius for different elapsed times, t e . In the lower 
panel the mass accretion rates at three different radii are displayed 
as a function of time. 

show the mass accretion rate, tin = 4irr 2 pv r , as a function 
of radius for different t e , which is the elapsed time from 
the beginning of collapse. It is seen that the rarefac- 
tion wave generated by the inflow at the inner bound- 
ary propagates outward, triggering the infall of matter 
at large radii. After t e ~ 300ms the accretion rates at 
r < 500km become independent of radius. This implies 
that the flows in this region can be approximated by 
steady accretions. The region is actually expanded out- 
ward gradually. The lower panel shows the accretion 
rates at three different radii as a function of time. As 
pointed out right now, they coincide with each other af- 
ter t e ~ 300ms. Before this time, on the other hand, 
the accretion rate is higher at smaller radii. There ap- 
pears a peak at t e ~ 180ms , which is rather insensi- 
tive to the radius. From a c omparison with realistic 
simulations (jJanka et al.ll2007t) . we find that this time 
roughly corresponds to the core bounce. We hence re- 
fer to as the time elapsed from this point as the post- 



bounce time hereafter, i.e, t 



pb 



t, 



t 



e(p) 



where t 



e(p) 



denotes as the time of peak accretion rate. In the same 
panel we also show the mass that has flown into the in- 
ner boundary by the given time, to which we refer as the 
proto-neutron star (PNS) mass (Mp^g). In the bottom 
panel we present the time evolution in the plain of mass 
accretion rate and PNS mass. 
The evolutions of the chemical composition together 
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Fig. 3. — The evolutions of the chemical composition as well as density, velocity and temperature. In the left column the density and 
velocity are displayed for three different times. In the right column, the abundances of some representative nuclei are shown with the 
temperature for the same three times. 



with the density, velocity and temperature are displayed 
in Fig. [3] In the left column the density and velocity are 
shown for three different post-bounce times (t p b)- As the 
time passes, the density at a fixed radius decreases mono- 
tonically whereas the inflow velocity gets larger. The 
outward propagation of the rarefaction wave is also rec- 
ognized in the figure, which is exactly how the implosion 
of envelope proceeds. In the right column, the chemical 
composition and temperature are presented for the same 
tpb. The temperature at a fixed radius is in general a 
decreasing function of time. It is observed that heavy 
elements are advected inward. In addition, the changes 



in composition by nuclear burnings are also taken into 
account in this figure. 

We employ these results not only at the next step in 
the preparation of initial conditions, which we will de- 
scribe in the next section, but also for the simulations of 
the post-relaunch evolutions, the results of which will be 
presented in § [3l 

2.4. Step 2: search of critical luminosities 

The aim of this step is to construct the critical steady 
accretion flows with a standing shock wave for the mass 
accretion rates obtained in step 1. It is important here 
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to define unambiguously the critical point for a given 
accretion rate, since in this paper it is not meant for 
the flow with the luminosity, above which n o steady ac- 
cretion flow exists iBurrows &: Goshvl (|1993l ) . Instead we 
define it to be the flow, in which the stalled shock wave 
is actually relaunched within a certain time. This is be- 
cause the shock revival normally occurs owing to hydro- 
dynamical instabilities even in ID before the luminosity, 
above which no steady accretion is possible, is reached. 
Hence we determine the critical point hydrodynamically 
by following the growths of the instabilities for initially 
spherically symmetric and steady accretion flows. 

We hence adopt a two-step procedure. In the first 
step, we construct a sequence of spherically symmetric 
and steady accretion flows with a standing shock wave 
for a given mass accretion rate. We solve Eqs. ([I])-© m 
ID, dropping the Eulerian time derivatives. At the shock 
wave, we impose the Rankinc-Hugoniot jump conditions. 
The nuclear reactions and weak interactions that are de- 
scribed in §2.21 are fully taken into account. The outer 
boundary of the computational domain is set to 500km 
and the values of various quantities are taken from the 
results of Step 1 at the times that correspond to the given 
mass accretion rates (see the bottom panel of Fig. [5]). At 
the inner boundary, which is placed according to Eq. (|8]l 
with the post-bounce time obtained from the mass ac- 
cretion rate, we impose the condition that the density be 
p = 10 u g/cm. We cover the computational region with 
300 grid points. As the neutrino luminosity is increased, 
the location of the standing shock wave is shifted out- 
wards and at some point the steady solution ceases to 
exist. As mentioned already, however, we do not need to 
search that point, since the shock revival occurs earlier 
owing to the hydrodynamical instabilities. We do need 
to identify this point, to which we refer the critical point 
in this paper, in the second step. 

The hydrodynamical simulations are performed both 
in ID and in 2D in the second step. As mentioned just 
now, these computations are used to judge whether the 
spherically symmetric, steady accretion flows, which are 
obtained in the first step, induce shock revival by the 
hydrodynamical instabilities. The nature of the instabil- 
ities are different between ID and 2D: in ID radial oscil- 
lations become over-stabilized at some luminosity, which 
is lower than the one, at which the steady flow ceases 
to exist (jOhnishi et al.ll2006tlFerni^d cz 2012); in 2D, on 
the other hand, the non-radial instability called S ASI oc- 
curs even earlier on (jYamasaki fc Yamada!l2007h . Hence 
in reality the latter will be more important. We think 
that ID models arc still useful to understand the physical 
processes that occur after the shock relaunch as well as 
to elucidate the differences caused by the dimensionality 
of hydrodynamics. 

We solve numerically Eqs. ([I]) - © with all the time 
derivatives turned on. Both the input physics and radial 
grid are identical to those employed in the first step, in 
which the steady accretion flows are calculated. In 2D 
simulations we deploy 60 grid points in the 6— direction 
to cover 180° and add random 1% perturbations to the 
radial velocity to induce SASI. In these simulations we 
fix both the outer and inner boundary conditions and fol- 
low the evolution for 200ms. If the shock wave reaches 
the outer boundary located at r = 500km within this 
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Fig. 4. — The critical neutrino luminosities in ID and 2D as a 
function of the post-bounce time. 

period, we judge that the shock is successfully revived. 
The reason why we fix the boundary conditions is that 
if the shock revival occurs at a certain time, the insta- 
bilities should have reached the nonlinear stage by that 
time but the growth of the instabilities takes some time. 
For each mass accretion rate we determine the minimum 
luminosity for the successful shock revival within a few 
percent and refer to it as the critical luminosity. In ID 
the shock relaunch is preceded by the growth of radial 
oscillations (not shown in the figure) whereas the non- 
radial modes with I = 1, 2 are responsible for the shock 
revival in 2D. Here I stands for the index of spherical har- 
monics in the expansion of unstable modes. In Fig. |4] the 
critical luminosities are presented both for ID and 2D as 
a function of mass accretion rate. It is evident that the 
critical luminosity is decreasing function of mass accre- 
tion rat e and it is reduced in 2D, both of which are w ell 
known (jOhnishi et al.ll2006l iMurphv fc Burrowsl[2008l ). 

2.5. Step 3: computations of post-relaunch evolutions 

In this section we give some details of the ID and 2D 
hydrodynamical simulations of post-shock-revival evolu- 
tions. We continue the computations of step 2 for the 
models with the critical luminosities. We first map the 
results to a larger mesh that covers the region extending 
from the neutrino sphere to the radius of r = 2 x 10 5 km. 
In all ID models we computed the post-revival evolutions 
for ~ 2s, which is found to be long enough to estimate 
the explosion energy. In fact, we follow the evolutions for 
two ID models until the shock reaches the stellar surface, 
which is located at r = 5 x 10 8 km. In those simulations 
we expand the mesh twice as the shock propagates out- 
ward. The inner boundary is also shifted to larger radii 
so that we could avoid too severe CFL conditions on the 
time step. In all simulations we employ non-uniform 650 
radial grid points in ID and 500 points in 2D. In 2D sim- 
ulations, 60 grid points are distributed uniformly in the 
direction to cover the entire meridian section. 

The outer boundary condition poses no problem this 
time, since it is located at a very large radius. We just 
impose the free in- flow/out- flow condition there. The 
inner boundary conditions are a bit more difficult. We 
assume the time evolution of neutrino luminosity is given 

by 

Lv{tex P ) = Lv,c ■ cxp(-i e;E p/800ms), (11) 

where t exp is the time elapsed from the shock relaunch 
and L„ jC is the critical luminosity obtained in Step 2. We 
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fix the density, pressure and velocity at the ghost mesh 
point at the inner boundary when matter is flowing in- 
ward. When matter begins to flow outward, i.e. the tran- 
sition to the neutrino wind phase occurs, those quantities 
are extrapolated from the innermost active mesh point 
to the ghost mesh point except when the entropy per 
baryon tends to be too high, in which case we put the 
upper bound of s = 100/cs to the entropy per baryon and 
the density is adjusted. These prescriptions are applied 
to each angular grid point at the inner boundary for 2D 
simulations. 

We investigate seven ID models, for which the stalled 
shock is relaunched at t pb = 200, 300, 400, 500, 600, 700 
and 800ms. Note that t p b has a one to one correspon- 
dence with the mass accretion rate, which is shown in 
Fig. [2J Five 2D simulations are performed, in which the 
shock revival is assumed to occur at t p b = 200, 300, 400, 
500 and 600ms. See Fig. [4] for the critical luminosities 
in these models. The input physics, such as nuclear and 
weak interactions, are the same as those employed in the 
second step of Step 2. The results of all the computa- 
tions in this step will be presented in the next section 
first for the 1 D models and then for the 2D cases. 

3. RESULTS 
3.1. Spherically symmetric ID models 
3.1.1. The evolution of the fiducial model 

We first describe in detail the evolution of the fD 
model, in which the shock relaunch is assumed to oc- 
cur at tpi, = 400ms. This corresponds to the time, at 
which the mass accretion rate is 0.52M Q /s. The critical 
luminosity in fD is 5 x 10 52 erg/s. 

In Fig. [5] we show the density, temperature and mass 
fractions of representative nuclei as a function of radius 
for four different times. In the upper left panel the pro- 
file at t exp = 20ms after the shock relaunch is displayed. 
The shock is still located around r = 500km. The post- 
shock temperature is T > lMeV, so high that the nuclei, 
mainly 28 Si, flowing into the shock are decomposed to a 
particles, which are further disintegrated into nucleons 
immediately. The post shock composition is perfectly 
described by NSE. At t exp = fOOms the shock reaches 
r ~ 2, 000km but is still inside the Si layer as seen in 
the upper right panel. The post-shock matter is mainly 
composed of a particles, which are not disintegrated any 
more owing to the lower temperature, T ~ 7 x 10 9 K. An- 
other 50ms later, the shock enters the Oxygen layer (see 
the lower left panel). Now 56 Ni emerges just behind the 
shock wave. This is mainly due to the recombination of a 
particles, which will be evident shortly. The post-shock 
temperature is T ~ 5 x 10 9 K and matter is beginning to 
be out of NSE. In the lower right panel we present the 
profile at t exp = 500ms. At this time, the temperature 
is T ~ 2 x 10 9 K and matter is completely out of NSE 
and the nuclear reactions yield mainly 28 Si. Much be- 
hind the shock wave some a particles are recombined to 
56 Ni. Slightly later all nuclear reactions are terminated 
behind the shock wave, since the temperature does not 
rise to high enough values by shock heating. 

In Fig. |6] we show the masses of proton, neutron, a, 
28 Si and 56 Ni integrated over the region inside the shock 
wave as a function of t exp for the ID fiducial model. 
In accord with the description in the previous para- 



graph, a particles are the main yield of nuclear reactions 
up to t exp ~ 100ms. The depletion of neutrons after 
texp ~ 50ms implies that the nucleons are recombined 
to a particles during this period. From t exp ~ f 00ms to 
texp ~ 150ms, on the other hand, a particles are dimin- 
ished while 56 Ni is increase, which means that the former 
is recombined to the latter. After t exp ~ f 50ms a parti- 
cles cease to recombine any more and are frozen, and 56 Ni 
and later 28 Si are produced by nuclear burnings. These 
results are obtained with the nuclear network with 28 nu- 
clei (see W2.2\i . In order to confirm that it is large enough, 
we conduct a larger network with 463 nuclei as a post- 
processing calculation, employing the time evolutions of 
density, temperature and electron fraction obtained by 
the simulation with the original network. The nickel and 
silicon masses are O.MOMq and 0.068M Q for the larger 
network, whereas they are 0.151M Q and O.O71M0 for the 
standard case. Both results are in good agreement, im- 
plying that the original network is appropriate. 

The energy release by these nuclear reactions is pre- 
sented in Fig. [7J We time-integrate the energy produc- 
tion rate over the entire evolution as a function of mass 
coordinates. In so doing, we distinguish the contributions 
of the recombinations from those of the nuclear burnings. 
Moreover, the former is divided into two pieces, one of 
which comes from the recombinations in NSE and the 
other from those in non-NSE. In the interior (< 1.75M Q ) 
the recombinations start in NSE and end in non-NSE and 
there are hence two contributions. In the slightly outer 
layer up to ~ 1.8M Q , on the other hand, the recombi- 
nations occur in non-NSE conditions. Further outside 
(> 1.8Mq) the nuclear burnings take their places. It is 
evident that the largest energy release comes from the 
recombinations that occur in NSE and the contributions 
of the nuclear burnings arc rather minor even after the 
integrated over the mass coordinate for this particular 
model. This is a generic trend as will be shown later. 
There is a gap between ~ 1.75M and ~ 1.8M Q . This is 
a region, in which energy is not released but absorbed. 
The main reaction in this region is the burning of 28 Si to 
56 Ni. Some fractions of 28 Si are disintegrated to a par- 
ticles, however. Although the latter is minor, the mass 
difference between a and 28 Si is much greater than that 
between 28 Si and 56 Ni. As a results the energy suck by 
the decomposition overwhelms the energy release by the 
burning. 

3.1.2. The evolution of diagnostic explosion energy 

Understanding the evolutions of density, temperature, 
chemical composition as well as the energy generations 
by nuclear reactions, we now turn our attention to the 
explosion energy. Following the conventional practice, 
we define the diagnostic explosion energy of provisional 
ejecta. At each grid point the total energy density, etot, 
is given by 

where ekin = l/2pv 2 is the kinetic energy density, e^t 
denotes the internal energy density and e grav = p($+& c ) 
stands for the gravitational potential energy density. We 
judge that the mass element at a certain grid point will 
be ejected if the total energy density is positive, e tot > 0, 
and if the radial velocity is positive (v r > 0) at a given 
time. Then the diagnostic explosion energy is defined 
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Fig. 5. — The density, temperature and mass fractions of representative nuclei for the ID fiducial model as a function of radius at four 
different post-relaunch times 
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Fig. 6. — The masses of proton, neutron, a, 28 Si and 56 Ni 
integrated over the post-shock region as a function of the elapsed 
time for the ID fiducial model. 



as a function of time to be the sum of the total energy 
density times volume over the ejecta, which is just deter- 
mined. 

The diagnostic explosion energy changes in time in- 
deed. In the early phase of shock revival, the neutrino 
heating is the main source of the diagnostic explosion 
energy. As the shock propagates outward, the neutrino 
heating becomes inefficient, since the matter to be heated 
is also shifted to larger radii, where the neutrino flux is 
lower, and the luminosity itself becomes smaller as the 



1.6 1.7 1.8 1.9 2 2.1 2.2 

M r [M ] 

Fig. 7. — Total energy release by nuclear reactions as a function 
of mass coordinates. The energy production rate is integrated over 
the entire evolution. Three contributions, i.e., recombinations in 
NSE (red line), those in non-NSE (left portion of green line) and 
nuclear burnings (right portion of green line), are distinguished. 



time passes. It is also important that nucleons, which 
are mainly responsible for the heating, are depleted as 
they recombine to a particles and heavier nuclei as the 
temperature decreases. After the neutrino heating sub- 
sides, the nuclear reactions are the main energy source. 
As described in the previous section, the recombination 
of nucleons occurs at first and the nuclear burnings take 
their place later. After all nuclear reactions are termi- 
nated owing to low temperatures, the diagnostic explo- 
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sion energy decreases slowly since matter, which is grav- 
itationally bound and hence has negative specific energy, 
is swallowed by the shock wave. As the shock wave pro- 
ceeds outwards, this contribution becomes smaller and 
the diagnostic explosion energy approaches its asymp- 
totic value, the actual explosion energy. 

le+51 | 1 1 , 1 1 




Fig. 8. — The time evolution of diagnostic explosion energy 
for the ID fiducial model. The individual contributions from the 
neutrino heating (green line) as well as nuclear reactions (blue line) 
are also shown. 

Figure [5] shows the time evolution of the diagnostic ex- 
plosion energy for the ID fiducial model, in which the 
shock is relaunched at t p b = 400ms. The horizontal axis 
in the figure is the time elapsed from the shock relaunch, 
t exp . The diagnostic explosion energy increases for the 
first ~ 200ms. Then it decreases gradually and becomes 
almost constant at t exp ~ Is. Also displayed in the figure 
are the individual contributions to the diagnostic explo- 
sion energy from the neutrino heating and nuclear reac- 
tions. As described in the previous paragraph, the neu- 
trino heating is dominant over the nuclear reactions ini- 
tially up to t exp ~ 120ms. Then the nuclear reactions be- 
come more important and raise the diagnostic explosion 
energy to ~ 10 51 erg at by the time t exp ~ 200ms in this 
particular case. As indicated by the colored shades in 
the figure, the nuclear reactions are mainly the recombi- 
nations until t exp ~ 150ms. The nuclear burnings follow 
until t exp ~ 300ms. The asymptotic value of the diag- 
nostic explosion energy is approached from above owing 
to the engulfing of matter with negative energy by the 
outgoing shock wave. 

In order to confirm that the final explosion energy has 
been already reached in the above computation, we con- 
tinue to evolve this model until the shock wave reaches 
the stellar surface. We shift both the outer and inner 
boundaries as mentioned in ij2.5l to avoid too severe CFL 
conditions at the innermost mesh point,. The result is 
presented in Fig. HI It is clear that the diagnostic ex- 
plosion energy is essentially constant for t exp > Is. Also 
shown in the figure is the result for another model, in 
which the shock relaunch is delayed until t p b = 800ms. 
The accretion late is ~ 0.23M Q /s and the critical lu- 
minosity is ~ 3.3 x 10 52 erg/s for this ID model. As is 
obvious from the figure, the asymptotic explosion energy 
is considerably smaller, ~ 1.1 x 10 50 erg, and we have to 
wait for ~ 2s before the diagnostic explosion energy is 
settled to the asymptotic value. This is a generic trend: 
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Fig. 9. — The long-term evolutions of diagnostic explosion ener- 
gies for the ID fiducial model as well as for the model of the latest 
shock relaunch. 

as the shock relaunch is delayed, it takes more time to 
reach the final explosion energy. 

3.1.3. Systematics 
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Fig. 10. — The asymptotic values of diagnostic explosion energy 
for all ID models. 

In this section, we look into the results of other models 
and see how generic what we have found so far for the 
fiducial model is. In Fig.[TU]we show the asymptotic val- 
ues of diagnostic explosion energy for different models as 
a function of the shock-relaunch time. We can clearly see 
that the explosion energy is a monotonically decreasing 
function of the shock-revival time. This is main due to 
the fact that the mass of accreting matter gets smaller 
as the time passes and is nothing unexpected (see e.g. 
(|Scheck et al.l l2006j)). It is stressed, however, this is the 
first clear demonstration of the fact with the nuclear re- 
actions and EOS being taken into account consistently 
in sufficiently long computations, in which the diagnos- 
tic explosion energy is confirmed to reach the asymptotic 
value. 

Also shown in the figure are the individual contribu- 
tions to the diagnostic explosion energy from the nuclear 
reactions and neutrino heating. Both of them also de- 
crease as the shock revival is delayed. It is found, how- 
ever, that the contribution of the nuclear reactions dimin- 
ishes more rapidly. This is simply due to the fact that 
the temperature rise by the shock passage is smaller in 
weaker explosions. Note that the explosion energy are 
smaller than the sum of the two contributions, since the 
accretion of gravitationally bound matter gives negative 
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contributions to the explosion energy as mentioned al- 
ready. It should be also emphasized that the recombi- 
nation energy is eventually originated from the neutrino 
heating because the recombinations are necessarily pre- 
ceded by the endothermic dissociations of heavy nuclei 
that exist prior to collapse and those consumed energies 
are replenished by neutrinos. The neutrino heating also 
plays a vital role to push the post-bounce configuration 
to the critical point and further heat up matter until they 
become gravitationally unbound in the earliest phase of 
shock revival. 
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Fig. 11. — The individual contributions from the nuclear recom- 
binations in NSE and those in non-NSE as well as nuclear burnings 
to the explosion energy. 

In Fig.[TT]we further divide the contribution of nuclear 
reactions into those from the recombinations in and out 
of NSE as well as from the nuclear burnings. Roughly 
speaking, the re-assemble of nucleons to a particles oc- 
curs in the recombination in NSE and the further recom- 
binations to heavier nuclei proceed in the environment 
out of NSE. As can be seen, all the contributions again 
decline as the shock relaunch is delayed. Regardless the 
recombination that occurs in NSE is the greatest con- 
tributor. The recombinations, both in and out of NSE, 
decline more rapidly than the nuclear burning and the 
latter contributes more than the recombination out of 
NSE for the model, in which the shock revival occurs at 
the latest time (t p b = 800ms) and the weakest explosion 
is obtained. The reason why the nuclear burning declines 
more slowly is that the temperatures obtained by shock 
heating is roughly proportional to the quarter power of 
the explosion energy. 

Next we turn our attention to the synthesis of 56 Ni, 
one of the most important observables in the supernova 
explosion. The synthesized mass of 56 N i is correlated 
with the explosion energy (jHamuvl 12003V. the greater 
the explosion energy is, the more 56 Ni is produced. This 
is demonstrated in Fig. [12] In the upper panel we again 
present the asymptotic values of the diagnostic explosion 
energy, which we simply refer to as the explosion energy 
here, as a function of the mass accretion rate at shock 
relaunch. The corresponding shock-relaunch times are 
given in the figure. In the lower panel, the mass of 56 Ni 
in the ejecta is displayed also as a function of the mass 
accretion rate at shock relaunch. The ejecta was defined 
earlier to be the collection of the mass elements that 
have positive total energy density and radial velocity (see 
fl3,1.2[) . The positive correlation between the explosion 
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Fig. 12. — The explosion energies and 56 Ni masses for all ID 
models. 

energy and the mass of 56 Ni in the ejecta are evident. 

In fact, 56 Ni may be produced too much in these ID 
models. The canonical explosion energy (~ 10 51 erg) is 
attained only by the models that relaunch the stalled 
shock wave relatively early (t p b < 400ms). On the other 
hand, the masses in the ejecta of 56 Ni synthesized for 
all these models are > 0.15M Q , which is substantially 
larger t han the values estim ated from observations, < 
O.1M (jSmartt et al.ll2009a[) . Note that the mass of 56 Ni 
ejected by SN1987A is estimated to be ~ 0.07M Q and 
corresponds to the shock-relaunch time of t p b ~ 600ms 
in our ID model; this rather late shock revival gives only 
a weak explosion of ~ 0.4 x 10 51 erg, smaller than the 
most likely ex plosion energy (0.9 x 10 51 erg) derived ob- 
servationally (Kas cn fc Wooslevl 120091 ). It is true that 
both of the observational estimates and the theoretical 
predictions presented here have uncertainties. In fact, 
the neutrino transport as well as evolutions of proto- 
neutron star, which are neglected and roughly mimicked 
in this paper, are the main source of uncertainties in 
the results shown above. We believe, however, that the 
general trends would be unchanged even if more sophis- 
ticated treatments were adopted. The above argument 
hence may be regarded as yet another reason that we do 
not believe that the ID neutrino heating works. 

To better understand the origin of the overproduction 
of 56 Ni we have done an ordinary calculation of explosive 
nucleosynthesis as a post-process for the densities, tem- 
peratures and electron fractions obtained for the mod- 
els presented above. In so doing, the so-called thermal 
bomb method, in which thermal energy is de posited ini- 
tially in the innermost region, is employed (jHashimotol 
Il995| ). The explosion energy and mass cut are chosen 
so that they agree with those of the original models. 
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TABLE 1 

Comparison with an ordinary explosive-nucleosynthesis calculation. 



Shock relaunch time Explosion energy Proto-neutron star mass 56 Ni mass 

[ms] [10 51 crg] [Mq] [M ] 

ID model thermal bomb ID thermal bomb ID thermal bomb 
300 E50 L50 L~60 L59 021 015 

400 0.95 .94 1.67 1.64 0.15 0.086 

500 0.50 0.60 1.70 1.73 0.10 0.043 




2 3 4 5 

Tpeak[10 9 K] 



Fig. 13. — Comparison of the distributions of peak temperature 
between the ID fiducial model and corresponding thermal-bomb 
model. 

Interestingly the calculation of explosive nucleosynthe- 
sis consistently produces smaller amounts of 56 Ni, which 
is given in Tabic 13.1.31 In fact, the fiducial model re- 
produces the observational estimate for SN1987A much 
better, which is just a coincidence though. Figure [13] 
compares the distributions of peak temperature between 
the ID fiducial model and corresponding thermal-bomb 
model. It is clear that the fiducial model has a larger 
amount of mass elements that achieve high enough tem- 
peratures to produce 56 Ni. It seems to take the neutrino 
heating mechanism a greater thermal energy to unbound 
the accreting envelope. This result may also be a caution 
in employing the thermal bomb method in the explosive 

nucleosynt hesis calculation s. 

In fact, lYoung k. Frverl ()2007l ) discussed uncertainty 
in nuclcosynthetic yields with ID spherically symmetric 
explosion models and showed that 56 Ni evaluated with 
the piston-driven model (E cxp = 1.2 xl0 51 erg) is larger 
by a factor of 2-3 than that with the thermal-bomb model 
(£cxp = 1.5 x 10 51 erg) with the same remnant mass. Our 
ID neutrino-driven explosion models may be closer to the 
piston-driven explosion model than to the thermal-bomb 
explosion model. 

3.2. Axisymmetric 2D models 

In the following we focus on the effect of dimension. In 
£12.41 we already showed that the critical neutrino lumi- 
nosity is lowered in the axisymmetric 2D models than in 
the spherically symmetric ID model for the same accre- 
tion rate at the shock relaunch. Performing 2D simula- 
tions further in Step 3, we are concerned with how and 
how much the results we obtained for ID models so far 
are modified in 2D cases. 

The input physics for 2D simulations is essentially 
the same as for the ID model except for the perturba- 
tions added to the radial velocities with random mag- 
nitudes up to 1% at the beginning of Step 2. We 



have investigated the models, in which the stalled shock 
wave is relaunched at the post-bounce times of t p b = 
200, 300, 400, 500 and 600ms. The computational domain 
is extended radially from r = 500km to r = 20, 000km in 
Step 3. 

3.2.1. Dynamics of aspherical shock revival 

We first look at the post-relaunch dynamics of the 
2D fiducial model, which is a 2D counterpart to the 
ID fiducial model, in which the stalled shock revives at 
tpb = 400ms. As shown in Fig. [TU which displays the 
contours of entropy per baryon and mass fraction of 56 Ni, 
the shock expansi on is highly asp h erical , which was also 
demon s trated in iKifonidis et al.l ()2006l ): lOhnishi et al.1 
(120061) : IScheck et al.l (120061 e.g.). The shock front is 
elongated in the direction of the symmetry axis. The 
shock propagates more rapidly in the northern hemi- 
sphere whereas matter attains higher entropy per baryon 
in the opposite hemisphere. These features conform with 
the dominance of £ = 1, 2 modes owing to SASI. Here £ 
stands for the index of Legendre polynomials, which are 
included in the eigenfunctions of linearly unstable modes. 

As mentioned in £12.41 the shock revival occurs at a 
lower luminosity in the 2D model (L c = 4.2 x 10 52 erg/s) 
than in the ID counterpart (L c = 5.0xl0 52 erg/s). Unlike 
the ID case, some matter continues to accrete, forming 
down drafts particularly in the equatorial region, until 
much late times after the shock revival. As a consequence 
the (baryonic) mass of neutron star is larger in the 2D 
case (~ 1.8M Q ) than in the ID case (~ 1.65M©). This 
is actually a generic trend as shown later (see Fig. H7| . 
Another interesting feature found in the 2D model is the 
distribution of maximum temperatures that each mass 
element attains, which is obtained from the Lagrangian 
evolutions of tracer particles distributed in the ejecta. 
Figure [TS] shows the result in a histogram. It is evi- 
dent that the mass that reaches T — 5 x 10 9 K is larger 
in ID than in 2D. In the same figure we also show the 
distribution of electron fraction, y e (NSE), which is esti- 
mated when T becomes the boundary value of NSE, or 
7 x 10 9 K. Note that F e (NSE) is useful for the diagnosis 
of nuclear yields in ejecta. In the ID case, there exist 
too massive ejecta with Y^(NSE) < 0.49, which will pro- 
duce unacceptable amount of neutron-rich Ni isotopes 
an d 64 Zn compared with the solar abundances as shown 
in iFuiimoto et ail (|2011l ). Their overproduction of the 
slightly neutron-rich ejecta with Y e (NSE) < 0.49 in the 
ID case disappears in the 2D case owing to more ef- 
ficient neutrino interactions. Multi-dimensional models 
are therefore preferable in the point of view of the Galac- 
tic chemical evolution of isotopes although they may de- 
pend on the treatment of neutrino transfer. 

The post-relaunch evolutions of other models are qual- 
itatively similar to that of the fiducial model described 
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Fig. 14. — Contours of entropy (left panel) and mass fraction of 56 Ni (right panel) at t exp = 400ms for the 2D model, in which the 
stalled shock revives at t v h = 400ms. 
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Fig. 15. — Comparison of the distributions of peak temperatures Tp Ca ]j (left panel) and y e (NSE) (right panel) between the ID fiducial 
model and the 2D counterpart. 



just now. SASI is dominated by the £ = 1,2 modes again, 
making the shock front generically prolate with a marked 
equatorial asymmetry, which are displayed in Fig. 1161 
The accretion continues until long after the shock is re- 
launched. The resultant mass of neutron star is larger 
in 2D than in ID as pointed out already for the fidu- 
cial model and now shown in Fig. [TTl As the shock re- 
vival is delayed and the critical luminosity gets lower with 
the mass accretion rate beings smaller, the post-relaunch 



evolution becomes slower as in ID. This is even more so 
for the 2D models, for which the critical luminosity is 
smaller than for ID owing to the hydrodynamical insta- 
bilities (see Fig. 2]). How these hydrodynamical features 
affect the explosion energy and nickel mass is our primary 
concern and will be addressed in the next section. 

3.2.2. Diagnostic explosion energies and masses of Ni in 

the ejecta 
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Fig. 16. — Post-shock-relaunch distributions of entropy per baryon in the meridian section for all the 2D models at t cxp ~ 0.25 s. 



The diagnostic explosion energy is shown in Figure [18] 
as a function of t exp for the 2D fiducial model, in which 
the stalled shock wave is relaunched at t p b = 400ms. Also 
presented are the individual contributions from the neu- 
trino heating and nuclear reactions. As in the ID fidu- 
cial model (see Fig.HJ, the neutrino heating is effectively 
closed at t exp ~ 100ms. By this time the shock front 
has reached the location of r > 2,000 km, which is far 
enough from the proto-neutron star for the neutrino flux 
to become negligibly small. After the freeze-out of neu- 
trino heating, the diagnostic explosion energy increases 
via the nuclear reactions such as the recombination of 
4 He to heavier nuclei in the early phase and the Si and 



O burnings later on as shown in Fig. 1191 in which the 
masses of n, p, 4 He, and a particles, 28 Si and 56 Ni are 
shown as a function of time. The diagnostic explosion 
energy becomes almost constant in t exp ~ 300ms in this 
model. 

These features are essentially the same as what we ob- 
served for the ID counterpart (Figs. [B]and[5} except for 
the fact that the proportion of individual contributions 
is different: the nuclear reactions are overwhelmed by 
the neutrino heating in this particular case. Also shown 
in the figure (thin solid line) is the time variation of di- 
agnostic explosion energy of the simulation, in which all 
nuclear reactions are manually turned off in the non-NSE 



14 



Yamamoto et al. 




fpb(s) fpb(s) 



Fig. 17. — Comparison of the masses of proto-neutron stars at 
the end of computations between the ID and 2D models. 




Fig. 18. — The time evolution of diagnostic explosion energy 
for the 2D fiducial model, in which the shock is relaunched at 



: 400ms. 
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Fig. 19. — The time evolution of masses of n, p, He, Si, and 
56 Ni as a function of t exp for the 2D fiducial model, in which the 



shock is relaunched at t 



pb 



400ms. The evolution is similar to 



that in the ID counterpart (Fig[6]l. 

regime, i.e., T < 6 x 10 9 K. This computation is meant 
to distinguish the contribution of the recombinations in 
NSE from those of the recombinations that occur in non- 
NSE as well as those of nuclear burnings. The difference 
in the asymptotic diagnostic explosion energy between 
the two cases is 0.17 x 10 51 erg, roughly half the total 
contribution of nuclear reactions, 0.36 x 10 51 erg. 

Figure [20] shows the explosion energies for all the 2D 
models (thick lines) in comparison with those for the ID 
models (thin lines). It is interesting that not only the 
general trend that the explosion energy is a monotoni- 



FlG. 20. — Comparison of the explosion energies between the 
ID and 2D models. 

cally decreasing function of the shock-relaunch time but 
the explosion energy for a given shock-relaunch time is 
also similar between the 2D and ID models. This implies 
that it is the mass accretion rate that primarily controls 
the explosion energy. It should be noted, however, that 
the critical luminosity for a give mass accretion rate is 
smaller in 2D than in ID (see Fig. [4]). For a given neu- 
trino luminosity, the explosion energy is hence larger in 
2D. This is the advantage of non-spherical hydrodynam- 
ics that is commonly mentioned in the literature. 
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Fig. 21. — The individual contributions of neutrino heating and 
nuclear reactions to the explosion energy for all the 2D models. 



In Fig. [5T] we present the individual contributions to 
the explosion energy from the neutrino heating and nu- 
clear reactions, which should be compared with Fig. 1101 
The most remarkable is that the contribution of the 
nuclear reactions is much reduced in 2D and is domi- 
nant over that of the neutrino heating only for the mod- 
els with an early shock-relaunch (t p t, < 300ms). Al- 
though the contribution from the neutrino heating is 
also smaller, the difference between ID and 2D mod- 
els is much smaller. The reduction of the contribution 
from the nuclear reactions is directly related with the de- 
crease of the mass elements in the ejecta that attain high 
peak temperatures, which we have pointed out already 
for the fiducial model (see Fig. [T5|) . These reductions 
arc compensated for by the reduction in the (negative) 
contribution from the gravitational energy of accreting 
matter that is swallowed in the ejecta. The latter re- 
duction is a consequence of the fact that the accretion 
continues in 2D even after the shock is relaunched and 
not all the accreting matter is added to the ejecta. 
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Fig. 22. — Comparison of 56 Ni masses in the ejecta between the 
ID and 2D models. 

The reduction of the contribution to the explosion en- 
ergy from the nuclear reactions is also reflected in the 
production of 56 Ni, which is demonstrated in Fig. [5U It 
is apparent that the mass of 56 Ni is always smaller in 
2D than in ID as a function of the shock-relaunch time. 
It should be reminded that the ID models tend to over- 
produce the nickels; as discussed in §3.1.3| if the typical 
mass of 56 Ni in the ejecta of C CSNe is M Ni 5 0AM Q as 
observations seem to indicate (|Smartt et al.ll2009al) . the 
ID models require that the shock should be relaunched 
later than t p b ~ 500ms ; this implies a rather weak explo- 
sion (E exp < 0.5 x 10 51 erg), however; no ID model hence 
can give both the explosion energy and nickel mass in the 
appropriate range. In the 2D models, on the other hand, 
this problem is much relaxed. Indeed the explosion en- 
ergy is large enough if the shock relaunch occurs earlier 
than t p b ~ 500ms and the nickel is not over-produced if 
the shock is revived later than t p b ~ 300ms. Although it 
is entirely a different issue whether and how the critical 
luminosity is obtained, this range of shock- relaunch time, 
tpb ~ 300 — 500ms, may be regarded as the appropriate 
time for shock revival in fact. It is nice that the 2D mod- 
els have an "allowed" range, since the hydrodynamics is 
inevitably non-spherical owing to the hydrodynamical in- 
stabilities. Whether three-dimensional hydrodynamics, 
which is the reality, alters the result for 2D will be an 
important issue and will be studied in the future. 

Comparison between the abundances in SN ejecta and 
those in the so lar system will possibl y lead to the similar 
allowed range. iFuiimoto et al.l (|2011h performed detailed 
nucleosynthctic calculations f or the ejecta of SN explo - 
sion from a 15M Q progenitor (|Wooslev fc Weaverlll995l) . 
They were based on simulations of neutrino-driven as- 
pherical explosion, which employed with a hydrodynamic 
code, that is similar to the one used in the present 
study but neglecting the energy generation through nu- 
clear reactions. They showed that the explosions with 
tpb ~ 200 - 300ms give £; eX p and M ( 56 Ni) in the allowed 
range and that the abundances in the ejecta are similar 
to those in the solar system. Detailed nucleosynthesis 
studies taking into account of the energy generation via 
nuclear reactions will be interesting, since the feedback, 
which elevates entropy in the ejecta, will possibly en- 
hance the amounts of 56 Ni a nd 44 Ti, which a re slig htly 
and highly underproduced in iFujimoto et al.1 ()2011l ). re- 
spectively. 

4. SUMMARY AND DISCUSSION 



We have investigated by numerical experiments done 
in ID and 2D the post-shock-relaunch evolutions in the 
neutrino heating mechanism. Taking into account the 
fact that the shock revival should occur somewhere on 
the critical line in the L v — M diagram but exactly 
where is rather uncertain theoretically at present, we 
have treated the luminosity (or cquivalently the mass 
accretion rate) at the shock relaunch as a free parame- 
ter that we vary rather arbitrarily. Only the post-bounce 
phase has been computed and we have discarded the cen- 
tral region (r < 50km) and replace it with the prescribed 
boundary conditions. We have also neglected the neu- 
trino transport entirely and employed the light bulb ap- 
proximation. The shock revival is hence induced by giv- 
ing a critical luminosity at the inner boundary by hand. 

The critical luminosity itself has been determined 
also by hydrodynamical simulations, since not the non- 
existence of steady state but the onset of hydrodynam- 
ical instabilities dictates the shock revival. The mass 
accretion rate as a function of time is obtained by the 
simulation of gravitational implosion of a massive star 
envelope. We have adopted a r ealisti c 15M progenitor 
provided bv IWooslev fc Hegerl (|2007t ). In these compu- 
tations, we have taken into account the nuclear reactions 
among 28 nuclei that include 14 a nuclei as well as their 
feedback to hydrodynamics consistently. As a result, we 
have confirmed that the critical luminosity is a monoton- 
ically decreasing function of the shock-relaunch time (or 
cquivalently a monotonically increasing function of the 
mass accretion rate) and that 2D dynamics reduces the 
critical luminosity compared with ID dynamics. This 
is due to the non-radial hydrodynamical instabilities and 
the resultant enhancement of neutrino heating in the for- 
mer case. 

After re-mapping, we have continued the simulations 
until long after shock revival. In fact, for some ID models 
we have followed the post-shock-relaunch evolutions up 
to the shock breakout of the stellar surface, confirming 
that the diagnostic explosion energy has approached the 
asymptotic value much earlier and the shorter computa- 
tion times for other models are sufficient indeed. We have 
employed the same set of input physics as in the simu- 
lations for the setup of initial conditions. Integrating 
the source terms in the equation of energy conservation, 
we have evaluated the individual contributions from the 
neutrino heating and nuclear reactions to the explosion 
energy. In so doing, we further divide the latter to the 
contributions from the recombinations in and out of NSE 
as well as from the nuclear burnings. The axisymmetric 
2D simulations have been performed to elucidate the ef- 
fect of multi-dimensionality on the outcome. 

What we have found in these model computations are 
summarized as follows: 

1. Immediately after shock relaunch the neutrino 
heating is the dominant source of the explosion 
energy but is terminated before long as the shock 
propagates outwards. Then the nuclear reactions 
take its place, with the recombinations of nucleons 
to a particles under the NSE condition occurring 
first. As the temperature decreases, the NSE be- 
comes no longer satisfied. The recombinations of 
a particles to heavier elements proceed mainly in 
this non-NSE circumstance. When the tempera- 
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ture lowers further, the nuclear burnings of silicons 
and oxygens take place in the matter that flows 
into the shock wave. Matter that flows into the 
shock wave contributes negatively to the explosion 
energy, since it is gravitational bound. 

2. The final explosion energy is a monotonically de- 
creasing function of the shock-relaunch time (or 
cquivalently an increasing function of the mass ac- 
cretion rate at shock relaunch) irrespective of the 
dimensionality of hydrodynamics. There is no big 
difference between ID and 2D as long as the mass 
accretion rate is the same at shock relaunch. This 
implies that the mass accretion rate is the primary 
factor to control the explosion energy. Since the 
critical neutrino luminosity for a given mass accre- 
tion rate is lower in 2D than in ID, the explosion 
energy for a given neutrino luminosity is larger in 
2D. 

3. As the shock relaunch is delayed, it takes longer 
the diagnostic explosion energy reaches the asymp- 
totic value. In our ID fiducial model, in which the 
stalled shock is revived at t p b = 400ms and we ob- 
tain the explosion energy of E exp ~ 10 51 erg, the 
diagnostic explosion energy attains the asymptotic 
value in t exp ~ 300ms whereas it takes ~ 2s for 
the model, in which the shock relaunch is assumed 
to occur at t p b = 800ms and the explosion energy 
is as low as ~ 10 50 crg. The similar trend is also 
observed in the 2D models. 

4. In ID the nuclear reactions always overwhelm the 
neutrino heating. The difference becomes smaller 
as shock relaunch is delayed. To the nuclear reac- 
tions the recombinations of nucleons to a particles 
that occur mainly under the NSE condition are the 
dominant contributor, followed by the recombina- 
tions of a particles to heavier elements in the non- 
NSE condition. The nuclear burnings provide the 
smallest contribution except for the weakest explo- 
sion, which obtains for the latest shock-relaunch. 
In that case, the nuclear burnings beat the recombi- 
nations in non-NSE. In 2D, on the other hand, the 
neutrino heating and nuclear reactions give compa- 
rable contributions to the explosion energy with the 
former being larger in stronger explosions and vice 
versa. Note, however, that the rather small contri- 
bution of neutrino heating is deceptive in the sense 
that it is the ultimate source of the energy obtained 
from the nuclear recombinations and that it is cru- 
cial to set the stage for shock revival, which is not 
accounted for in the diagnostic explosion energy. 

5. In the ID models nickels are overproduced owing 
mainly to a larges mass that achieves high peak 
temperatures compared with the ordinary calcu- 
lation of explosive nucleosynthesis in post-process. 
In fact there is no ID model that gives the typi- 
cal values to the explosion energy and nickel mass 
simultaneously. Given observational and theoret- 
ical uncertainties, we are not certain how serious 



a problem this is for the moment. One may con- 
sider, however, that this is yet another reason to 
abandon the ID neutrino heating model. In 2D, 
on the other hand, this problem is solved, opening 
up the allowed region in the shock-relaunch time 
around t p b ~ 400ms , which produce the explosion 
energy and nickel mass in the appropriate range. 
This happens mainly because the mass of matter 
in the ejecta that attain high enough peak temper- 
atures is smaller in 2D. This is in turn related with 
the fact that the expansion and accretion occur si- 
multaneously in 2D, which is indeed reflected in the 
mass of proto-neutron star, which is larger in 2D 
at any post-bounce time. 

In the present paper we have employed the single 15Mq 
progenitor model, which we think is one of the most 
representative to produce the typi cal Type IIP CCSN. 
Very recently lUgliano et al.l (|2012l ) reported a possible 
stochastic nature in the outcome of the shock revival in 
the neutrino heating mechanism based on systematic ID 
hydrodynamical simulations. Although the stochasticity 
is less remarkable in the low mass end, it is hence manda- 
tory to extend the current work to other progenitors and 
see how generic our fin dings obtained in this paper are 
(|Yamamoto et al]|2012[ ). 3D models are also the top pri- 
ority in the future work, since we know that 3D SASI 
is qualitativel y different from 2D SASI we have studied 
in this paperj Iwakami et al.ll2008f) . It should be also re- 
called that lNordhaus et al.l (|20 iW claimed that shock re- 
vival is even easier in 3 D than in 2D althou gh controver- 
sies are still continuing (jHanke et al.ll2012h . We are par- 
ticularly concerned about how the allowed region in the 
shock-relaunch time that is opened in 2D is modified in 
3D. The relative importance of the nuclear reactions for 
the explosion energy compared with the neutrino heating 
is the highest in ID. We are certainly interested in what 
about 3D. 

One of the greatest uncertainties in the present study is 
the effect of the inner boundary condition that is imposed 
by hand. The eventual answer should come from fully 
consistent simulations of the entire core. The neutrino 
transport, which is neglected completely in this paper, 
is another factor that may modify the results systemati- 
cally. These caveats notwithstanding we believe that the 
results obtained in this paper are useful to understand 
the post-shock-relaunch evolution in the neutrino heat- 
ing mechanism, particularly how the diagnostic explo- 
sion energy approaches the final value. One of the goals 
of our project is to seek, probably after more system- 
atic investigations suggested above, the way to estimate 
the explosion energy from the early stage of post-shock- 
revival evolution, since realistic simulations may not be 
affordable for a few seconds after shock relaunch. 
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